RELATIVISTIC BURGERS EQUATIONS ON CURVED SPACETIMES. 
DERIVATION AND FINITE VOLUME APPROXIMATION 

PHILIPPE G. LEFLOCH*, HASAN MAKHLOF* , AND BAVER OKUTMUSTUPJ 

Abstract. Within the class of nonlinear hyperbolic balance laws posed on a curved spacetime (endowed with 
a volume form), we identify a hyperbolic balance law that enjoys the same Lorentz invariance property as the 
one satisfied by the Euler equations of relativistic compressible fluids. This model is unique up to normalization 
and converges to the standard inviscid Burgers equation in the limit of infinite light speed. Furthermore, from 
the Euler system of relativistic compressible flows on a curved background, we derive both the standard inviscid 
Burgers equation and our relativistic generalizations. The proposed models are referred to as relativistic Burgers 
equations on curved spacetimes and provide us with simple models on which numerical methods can be developed 
and analyzed. Next, we introduce a finite volume scheme for the approximation of discontinuous solutions to 
these relativistic Burgers equations. Our scheme is formulated geometrically and is consistent with the natural 
divergence form of the balance laws under consideration. It applies to weak solutions containing shock waves and, 
most importantly, is well balanced in the sense that it preserves steady solutions. Numerical experiments are 
presented which demonstrate the convergence of the proposed finite volume scheme and its relevance for computing 
entropy solutions on a curved background. 



1. Introduction. In computational fluid dynamics, the (inviscid) Burgers equation 

d t u + d x (u 2 /2) = 0, u = u(t,x), t > 0, x e M, 

is an important toy model, which can be formally derived from the Euler equations of compressible 
fluids 

d t p + d x (pu) = 0, d t {pu) + d x (pu 2 + p{p)) = 0, 

(p, u denoting the density and velocity of the fluid and p = p(p) its pressure) . The Burgers equation 
is also the simplest example of nonlinear hyperbolic conservation law and has played an essential 
role in the development of robust and accurate, shock-capturing schemes for the approximation of 
entropy solutions to nonlinear hyperbolic systems. 

In the present paper, we derive several relativistic generalizations of Burgers equation which, in 
addition, we formulate on a curved background, and we study the discretization of these equations 
via a well balanced, finite volume strategy based on time-independent solutions to these relativistic 
Burgers equations. More precisely, we consider the following class of nonlinear hyperbolic balance 
laws 

(1.1) div u (T(v)) = S{v), 

posed on an (n + l)-dimensional curved spacetime (with boundary), denoted by (M,u>) and en- 
dowed with a volume form ui. The unknown function is a scalar field v : M — > R and div^ denotes 
the divergence operator associated with the volume form. The flux vector field T — T(v) is a field 
of tangent vectors prescribed on M and depending smoothly upon v, while the right-hand side 
S = S(v) is a scalar field prescribed on M which also smoothly depends upon v. 

Hyperbolic conservation laws on curved spacetimes have been analyzed in recent years by 
LeFloch and collaborators; cf. the review [5], as well as [TJ[H[7]. In order to formulate the initial 



value problem for (1.1 1, we assume that the spacetime is foliated by hypersurfaces, that is, 



(1.2) M = |J H t , 



t>o 



in which each slice H t is an ?i-dimensional manifold, canonically endowed with a normal 1-form 
field N t and with the same topology as the one of the initial slice Hq. The global hyperbolicity 
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of the spacetime and the equation (1.1 1, by definition, requires that each slice H t be spacelike, in 
the sense that the function 



(1.3) v i — ^ T°(v) := (N t ,T(v)) is strictly increasing:. 



The class of nonlinear hyperbolic equations ( |l.l[ )-( 1.3 ) provides us with a scalar model on which 



one can develop and analyze numerical methods of approximations which are then applicable to 
the Euler equations of relativistic compressible fluid flows on a curved spacetime. We refer the 
reader to the follow-up work by LeFloch and Makhlof 6] which addresses this generalization. For 



the well posedness theory of the initial value problem associated with ( 1.1 ), wc refer the reader to 
P] and the references cited therein. 

In this paper, we introduce a geometric formulation of the finite volume method, which is 
classically formulated without taking the underlying geometry into account, and we formulate a 
well balanced version of this method, which is inspired from the work of Russo and collaborators 
[TU1 fTTl fT2l H"3"] and is built upon the Nessyahu-Tadmor scheme [9 . For a review on well balanced 
techniques, we refer the reader to Bouchut [3], Greenberg and Leroux [4], and LeVeque [8]. 

An outline of the present paper is as follows. In Section 2, we derive our relativistic version of 
Burgers equation when all geometric effects are neglected. In Section 3, we explain how to take the 
geometry into account and we also connect our model to the Euler equations. Next, in Section 4, 
we present the finite volume strategy, and finally, in Section 5, several numerical experiments. 
Section 6 contains concluding remarks on this paper. 

2. A Lorentz invariant conservation law. 

2.1. Derivation of the new model. In this section, we search for flux vector fields T = T(v) 



for which solutions to the equation ( 1.1 ) satisfy the same Lorentz invariant property as the solutions 
to the Euler equations of relativistic fluids. For the sake of simplicity in the presentation (and 
without loss of generality for the purpose of this section), we assume that n — 1, S(v) = 0, and 
that the manifold M = [0, +oo) x E is covered by a (single) coordinate chart (x°, x 1 ) in which the 



volume form reads w = dx°dx 1 . Hence, ( 1.1 ) takes the simplest form of a hyperbolic conservation 
law, i.e. 

d o T°(v)+d 1 T 1 (v) = 0, 

with do — d/dx° and di — d/dx 1 , while x° £ [0, +oo) and x 1 e K. In addition, in the present 
section, we also assume that the functions T° = T°(y) and T 1 = T 1 ^) are independent of 
(a; , a; 1 ). A formulation taking geometric effects into account will be introduced and investigated 
in Section |3j below. 

The Lorentz transformation (a; , a; 1 ) i— > (a; ,^ 1 ) reads 

(2.1) x° :=7 e (V) (x°-e 2 Vx 1 ), x 1 := j e (V) ( - V x° + a 1 ) , 

where e 6 (—1,1) denotes the inverse of the (normalized) speed of light, and 7 e (V) = (l — 

— 1/2 

e 2 V 2 ) the so-called Lorentz factor associated with a given speed V € (— 1/e, 1/e). These 
equations can be inverted to give 

x ° = 7e (y) (x° + eW), x 1 = le (V) (VxV + x 1 ). 

Instead of V, it is often convenient to use the normalized speed U £ K defined by 

, 1/2 

e eC/ := % (V)(l + eV) 




so that the Lorentz transformation takes the compact form x° ± ex 1 = e^ eU (x° ± ex 1 ) or, equiv- 
alently, 

(2.2) x° = cosh(eU)x° - e sinh(etf) a 1 , x 1 = -- sinh(eLT) x° + cosh(eC/) x 1 , 
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with cosh(e£7) = 7 e (V) and sinh(eJJ) = eV j e (V). Recall also that the set of Lorentz transforma- 
tions (together with the spatial rotations if n > 2) forms the so-called Lorentz isometry group, 
characterized by the condition that the length element of the Minkowski metric is preserved, that 

e~ 2 (x ) 2 + (a; 1 ) 2 . The relativistic Euler equations of compressible flu- 



1\2 



is, -e~ 2 (a; ) 2 + (x 1 ) 

ids are invariant under the Lorentz transformation (2.1). More precisely, given a speed V, the 
fluid velocity component v in the coordinate system (x u , x 1 ) is related to the component v in the 
coordinates (x^x 1 ) by the relation 



(2.3) 



V 



1 - e 2 Vv' 

Clearly, in the nonrelativistic limit e — > 0, one recovers the Galilean transformation 



(2.4) 



x 1 = -Vx° 



v = v — V when e = 0, 



in which the speed V belongs to R. 

Theorem 2.1 (Relativistic version of Burgers equation). The conservation law 



(2.5) 



d T°(v) + d 1 T 1 (v) = 



is invariant under the Lorentz transformations (2.2)-(2.3) if and only if, after suitable normaliza- 
tion (discussed below), one has 



(2.6) 



T» = 



vT 



t 2 v 2 ' 



- 1 



1 



where the unknown scalar field v takes its value in the interval (— 1/e, 1/e). 

We refer to (2.5)-(2.6) as the relativistic Burgers equations (on Minkowski spacetime). 

Proof. We first change coordinates and construct certain auxilliary functions, denoted below 
by H±(v). Precisely, in the coordinates (x° , x 1 ), we can differentiate the Lorentz expressions (2.1 1 
and obtain 



dx° 
dx 1 

and, by the chain rule, 
d T -- 
diT 1 = 



1 



Vl 



-e 2 V 2 



le(V), 



dx° 
dx 1 
dx 1 
dx 1 



-e 2 V le (V), 

7e(V) 



dT°dx° i dT° dx 1 
dx° dx° 
dT 1 dx° 



dx 1 dx° 
dT 1 dx 1 
dx° dx 1 dx 1 dx 1 



dT 1 



dT 1 

le(V)V + W 1,(V) 



By substituting these formulas in (2.5), the conservation law reads 



= d T° + d 1 T l = d { le {V) T° - e 2 7 e (V) VT 1 ) + d^-V le (V) T° + j e (V) T 1 ), 
which has the same structure as (2.5| in the coordinates (x^x 1 ) provided 

v-V 



(2.7) 
(2.8) 



7e (V) (T°(v)-e 2 VT 1 (v))=T°l Y 



le(V) ( - VT°(v) + T 1 ^)) = T 1 1 J- 




for some constants C\,C2- This leads us to the desired invariance property and it remains to 
determine the general expressions of T° jT 1 . 



We can always normalize the flux terms so that T (0) = T (0) = 0, and then observe that 
the constants C\ = —T°(—V) and C2 = — T l {— V) are determined as functions of V. Next, 



multiplying the second equation in (2.7) by ±e and summing up with the first one, we find 



7e (V) (1 T eF)(T° ± eT 1 )^) = (T° ± eT 1 ) ^ - (T° i^ 1 )^). 



By setting (T° ieT 1 ) =: these two equations read 
(2.9) H ± (v) -- 



l_±eV ( v-V 



l-e 2 vV 



H±(-V) 



Next, we define the new variables 



1 . /1 + ev 

u := — In 

2e \l-ev 



e (u):=u, e (tr):=V, 



so that 



If 



2eu 



<l>e{u), 



1 e 2«(u-U) _ 1 



e e 



2e(u-£/) + J 



= ^(u-c0» 



l±eV 
1 Te^ 



By inserting these expressions in (|2.9|), it follows that 
(2.10) 



H 



(«)) - ff±(0 e (O)) = e ±eC/ (i/±(0 £ ( M - CO) - ff±(</»e(-C/))). 



Abusing notation, we write H,cj),e instead of H±,<f> e ,e , respectively. We differentiate 



(2.10) with respect to u and obtain H' (</>(«)) (j>'(u) = e eU H' (4>(u — U))(j)'(u — U) and, next, with 
respect to U 

H"((f>(u - Uj) (4>'{u - U)) 2 = eH'(4>{u - U))<f>'(u - U) - 4>" {u - U)H' (<p(u - U)). 

By letting U = 0, it follows that H"(<j){u)) ((f>'(u)) 2 = eH' {<p(u))<j>' (u) - cf>" (u)H' (<j>(u)) . The 
substitution K(u) := H(<p(u)) allows us to rewrite the latter equation as K"{u) = eK'(u), which 
has the following general solution (in which c, c are constants): 

e eu 

K(u) = c h c. 



Alternatively, this latter relation can be derived by dividing (2.10) by u ^ and letting u — > 



Returning to H±, we deduce that there exist constants so that H±{<f>(u)) — ±e 1 (e ±tu — l) 



Recalling that H±=T°± eT 1 with H± (</»(«)) = H±(v), we find T° = 
from which we deduce 



and T 1 = 



H+-H- 
2e 



T°(v) = T°(cj ) (u)) = 
T 1 ( U )=T 1 (^)) = 



2e 



2e 2 



= - sinh(ew) = u + 0(e 2 u 3 ), 



:(a»h(eu)-l) = ^-+0(eV). 
2 



Observe that T° and T 1 are linear and quadratic, respectively, and that, in the limiting case e — > 0, 
we recover the inviscid Burgers equation. Moreover, substituting u — ~ In ( in T°f^(ti)j 

and T 1 (0(u)), we thus arrive at (2.6). □ 
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2.2. Hyperbolicity and convexity properties. The proposed relativistic Burgers equa- 
tion retains several key features of the relativistic Euler equations: 

• Like the Euler system, it has a conservative form. 

• Like the velocity component in the Euler system, our unknown v is constrained to lie in 
the interval (— 1/e, 1/e) limited by the (inverse of the) light speed parameter. 

• Like for the Euler system, by sending the light speed to infinity one recovers the classical 
(nonrelativistic) model. 

Theorem 2.2 (Properties of the relativistic Burgers equation). 



1. The map w — T°(v) = 



is increasing and one-to-one from (—1/e, 1/e) onto 



2. In terms of the new unknown wet, the equation {2.5) is equivalent to 

d o w + d 1 Mw)=0, f e (w) = ^(-1 + Vl + e 2 ^ 2 )- 



(2.11) 



The flux f e is strictly convex and, therefore, the conservation law (2.5) is genuinely non- 
linear in the sense that q"^o[^ — d w f e {w) is strictly increasing or strictly decreasing in 
the variable T°(v). 

3. In the nonrelativistic limit e — > 0, one recovers the inviscid Burgers equation 



(2.12) 

Proof. We set w := T°(v) 



d v T°(v) 



<9 w + <9i(u 2 /2) = 0, u = w^V 1 ) e 



VT- 



(1 - e 2 v 2 ) 3 / 



and take its derivative with respect to v 
>0, v e (-1/e, 1/e), 



so that T°(v) is increasing and one-to-one from (—1/e, 1/e) onto M. We get v — ^ 14 ^ 2 „,2 (after 
choosing the positive branch, for the sake of definiteness) and, substituting in (2.5), we arrive at 
(2.11 ). The flux is strictly convex, since 



/» = 



Furthermore, we have 



Vl + e 2 w 2 '' 

d v THv) = 
d v T°(v) 



= (l + eW)3/ 2 > °- 



yl + e 2 w 2 

which is positive for w > and negative for w < 0. Finally, we have 



1 



1 



T°(v) = - sinh(eu) = u + (9(e 2 u 3 ), T\v) = ^ ( cosh(eu) - l) = — + 0(e 2 u 4 ) 



When e approaches 0, one recovers the standard inviscid Burgers equation (2.12). □ 

2.3. The nonrelativistic limit. The Galilean transformation (x ,^ 1 ) i-> (a; , a; 1 ) is defined 
by (2.4), in which the speed V G R is a given parameter. The velocity component v in the 
coordinate system (a; ,^ 1 ) is related to the velocity v in the coordinates (a; ,^; 1 ) by the relation 
v = v — V . Recall that one recovers the Galilean transformations from the Lorentz transformations 
in the limit e — > 0. 

Theorem 2.3 (A derivation of the (nonrelativistic) Burgers equation). The conservation law 



(2.13) 



d o T°(v)+d 1 T 1 (v)=0 



is invariant under the Galilean transformations (2.3)— (2.4) if and only if the flux T° and T 1 are 
linear and quadratic functions, respectively. After a suitable normalization, one has T°(v) = v 
and T x {v) = v 2 /2. 



Proof. Using the chain rule, we find 

gT0 _dT^dx°^ dT^dx 1 _ dT° &T° 
dx° dx° dx 1 dx° dx° dx 1 



dT 1 

dx 1 



By substituting this in (2.13), it follows that 



d T°{v + V) + Bi (t 1 {v + V) - VT°(v + V)) = 0. 

Now, this equation has the same algebraic structure as ( |2.13 1, provided 

(2.14) T°(v + V) = T°(v) + C 1 , T\v + V)-VT°(v + V) =T 1 {v) + C 2 , 

where C\ and C% are constants depending upon V. 

Then, we can derive the first equation in (2.14) with respect to v, that is, T°(tJ+ V") = T®(v), 
from which we conclude that T® is periodic of period V. As a result, T° is a linear function of the 
form T°(v) — Cv + C for some constants C and C. Next, multiplying the first equation in (2.14) 
by V and adding up to the second one, we get T x (v + V) = T x (iJ) + V T°(v) + VC X + C 2 - Since 
T° is a linear function and deriving twice the given equation, it follows that 

T^v + V) -T\v) = CVv + VC + Vd +C 2 , 
TUv + V)-TUv)-CV = 0, 

tUv + v)-tUv) = o. 

Therefore, T~ is periodic with period V, which implies that the first-order derivative is linear 
and that T 1 is quadratic. After normalization, one obtains exactly Burgers equation, corresponding 
T°(v) = v and T 1 (w) = v 2 /2. □ 

2.4. Derivation from the relativistic Euler equations. 



Constant density. We present another derivation of (2.11), now from the relativistic Euler 
equations which read 



(2.15) 



d 



P(p) + P c2 v 2 



d [ (p + pc 2 ) 2 _ 



+ Pj +9iUp(p) + pc 2 )^ 

( v 2 
fdi ( P (p) + pc 2 )- 



p{p) 



o. 

= 0, 



where p and u denote the mass-energy density and the (suitably normalized) velocity of the fluid, 
respectively, while the pressure p = p(p) is given by an equation of state depending on the fluid 
and the constant c > represents the light speed. 

By formally setting the density p (and thus the pressure p(p)) to be a constant in the second 
Euler equation above, we obtain 



do 



•ft 



0. 



Introducing the change of unknown z = 1 _" 2t)2 or v = 1± ^ e 1 2 + z 4 £ ^-, with c = 1/e, we find 



d z + ^<9i 
2 e z 



1± \/l + 4e 2 z 2 ) = 0, 



which is precisely (2.11 ) if we choose the positive branch and suitably normalize the equation. 

6 



2.5. Derivation from a general covariant equation on Minkowski spacetime. We 

now adopt a more geometric standpoint, but yet we begin by the case that the metric is flat. 
Following the notation in the introduction, we consider a general conservation law posed on a 
spacetime (M,g) endowed with a Lorentzian metric g, i.e. 

(2.16) V a T a (v)=0 inM, 

where v ; M — > R is the unknown scalar field, and T a = T a (v) is a prescribed vector field 
depending on v as a parameter. By analogy with the Euler equations, we impose that, for every 
constant v, T a (v) is a unit vector field, that is, g a pT a {v)T^ (v) = 1. 

Specifically, we now assume that (M,g) is the (1 + l)-Minkowski spacetime, described in 
coordinates (a; , a; 1 ) = (t,x), by g = -(T 2 dt 2 + dx 2 . We obtain -e^ 2 T°(v) 2 + T 1 ^) 2 = 1, thus 
T 1 (w) = ±yl + e~ 2 T°(v) 2 and, after choosing the positive branch, 

d T° +dx(C+ v/l + e- 2 (T°) 2 ) = 0. 

Here, the constant C was added for convenience. Finally, introducing the change of variable 
w = e~ 2 T° and normalizing the constant to be C — —1, we again recover the proposed relativistic 



Burgers equation (2.11) 



3. Burgers equations with geometric effects. 

3.1. Relativistic Euler equations on Schwarzschild spacetime. 

Vanishing pressure on flat spacetime. Our standpoint is now different and we start 
directly from the Euler equations and formally assume that, in the relativistic Euler equations 



(2.15), the pressure vanishes identically, that is, 




do A +ft A =0, 



These two equations can be rewritten as 

(c 2 - v 2 )(d p + vdip) + P (2vd v + (v 2 + c 2 )d lV ) = 0, 
v(c 2 - v 2 )(d p + vdip) + p((v 2 + c 2 )d v + 2vc 2 d 1 v) = 0, 

and we recover the classical Burgers equation (in flat spacetime) 

(3.1) 9o« + 9i(« 2 /2) = 0, 

which, therefore, appears to be also a relevant model, even in the relativistic setting. 



Schwarzschild geometry. We now want to take geometric effects into account in (3.1 ). The 
scalar equation of interest is still sought in the form 

(3.2) V Q T Q (w)=0, v.M^R 

for some prescribed vector field T a (v). To determine the model of interest, we start from the 
Euler equations posed on a curved spacetime and, for definiteness, assume that the background 
manifold (M, g) is Schwarzschild spacetime. The latter geometry describes a spherically symmet- 
ric black hole solution of the vacuum Einstein equations given, in suitably chosen coordinates 
(a; , x , x 2 , x 3 ) = (ct,r,6,(p) adapted to the symmetry, by 

(3.3) g = - (l - ^\ <? dt 2 + (l - ^ 1 dr 2 + r 2 (d6 2 + sin 2 6 dtp 2 ) 



for t > and r > 2m where the coefficient to > is referred to as the mass. This spacetime is 
spherically symmetric, that is invariant under the group of rotations operating on the space-like 
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2-spheres given by keeping t and r constant. It is static, since the vector field dt is a time-like 
Killing vector, and is asymptotic to (flat) Minkowski in the limit r —> oo. In (3.3), there is 
an apparent (not physical) singularity at r — 2m which, however, is removed in the so-called 
Eddington-Finkelstein coordinates. 
The Euler equations read 

(3.4) V a T$ = 0, Tg = (pc 2 + P )u a uP+p g? , 

where Tp is the energy-momentum tensor of perfect fluids, p is the mass-energy density of the 
fluid, and u = (u a ) denotes its unit velocity vector, so g a p u a u@ = — 1. This model is supplemented 
by an equation of state for the pressure p = p{p) . We assume that solutions to the Euler equations 
depend only on the time variable t and radial variable r, and that the non- radial components of 
the velocity vanish, that is, (u a ) — (u (t, r), u (t, r), 0, 0). Since u is unit vector, we obtain 

(, 2 ™v )2+ (- 1 ) 2 



r I ' 1 — 2m 



It is convenient to introduce the velocity component v := c(l — — ) 1 ^y, so that 
We obtain 



rpOl _ fOl rpll _ (-1 _ '^ rr ^\ mil rp22 P(P) 



rpOCl _ ^ JlUl _ rplji rpll ( \ — — '— ) J 111 r £ Z ' 1 — 

' ( -i 2m \ 1 ' V ~ I ' 

with 

r"" :=-! r\ 1-":=— T 



u;, ._ C 2 p + p(p)v 2 / C 2 2 ~ Q1 C 2 p+_p(p) ~ n V 2 p+ P (p) 2 



Vanishing pressure on Schwarzschild spacetime. We suppose first that the pressure p 
vanishes identically, and the Euler system in 1 + 1 dimensions on a Schwarzschild background 
takes the simplified form: 



& f T 01 J + 9 r ((r - 2m) 2 f ") - 3m f " + m t-^f» 

or, equivalently, 

„ / r 2 p \ ,/ 2m\„ / flo \ 2r — 2m 

(3 6) 9 | ^ ) + fl ^ TON )i9 | 1,2/9 | + ~ + ^ y2r 

\ c 2 — v 2 J \ r J \ c 2 — v 2 J r 2 (c 2 — v 2 ) 

By combining these two equations together, we arrive at 

/ 2m\ „ m . „ 

d t v + (l - — J vd r v = (v 2 - c 2 ) 

or, equivalently, 



2 

(3.7) d t {r 2 v) + d r (r(r-2m) V -^ 



2 2 

= rv — mc . 



This provides us with a variant of Burgers equation which incorporates two features of Schwarz- 
schild spacetime, i.e. the effect of spherical symmetry and the mass of the black hole. We refer to 
it as Burgers equation on Schwarzschild spacetime. 
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3.2. Static solutions of Burgers equation on Schwarzschild spacetime. 

Static solutions. Static solutions will play a central role in order to design our numerical 



method. We search for i-independent solutions to the equation (3.7| 
(3.8) \d r (r{r -2m)« 2 j 



2 2 

rv — mc . 



Using the change of variable X := c — v and Y := (r — m)/m, we find ^S- = —2 x _ Y i and, after 
integration, X = K 2 yqjy, where K € (0, c) is an arbitrary constant. We can thus describe all 
steady solutions v — v s (r) by 

(3.9) v s (r) = ±Jc?-K*(l-—), K e (0,c) 



r / 



or, equivalently, 



c 2 



(3.10) Y^t = K ' ^G(0,c). 

r 

Conservative form. In fact, the above calculations suggest to check whether the equation 



(3.7) admits the conservative form 

d (av) + ^d y (a (« 2 -c 2 )) =0, 

and we find that y = y(r) and a = a(y) must satisfy 

a'(y) —2m dr 
ay - 



a{y) r 2 1 



Ira 



Solving these equations yields a — r J r 2m an d y = r + 2m In (r — 2m), and we obtain the following 
conservative form of Burgers equation on Schwarzschild spacetime 



< 3 -"> M^«UiM^<» 2 -->) = <>• 

Interestingly the variable y above is essentially the Eddington-Finkelstein coordinate. 

3.3. Relativistic Burgers equation with geometric effects. We now return to the ap- 
proach adopted in Section 2. We assume that the manifold is described by a single chart and, 
after identification, we set M = M. + x R. In coordinates (x°, a; 1 ) with d a := d/dx a (with a — 0, 1), 
the hyperbolic balance law under consideration reads 

(3.12) d Q (uJT {v, x°, x 1 )) +d 1 (oJT 1 {v, x ^ 1 )) ^uSiv^ ^ 1 ), 

where v : M — > K is the unknown function and T a — T a (v,x° ,a; 1 ) and S = S(v, x° , x l ) are 
prescribed (flux and source) fields on M, while uJ — ^(a; ,^ 1 ) is a positive weight-function. This 
equation is hyperbolic in the direction d/dx 1 provided d v T°(v) > 0, which is always assumed. 

For instance, in the case that the flux and source are independent of x°, x , the above balance 
law can be rewritten in the form 

d Q v + d 1 f(v) = s(v), d v f(v)- ,)rV{ " ] 



d v T°(v) : 



S(v) := ( S ( v ) - 9onT°(v) - d 1 nT\v)), 



:= In uj. 



It may occur that S vanishes identically, however, one should keep in mind that (3.12) is the 
geometric form of this equation: it should preferred for the numerical discretization. It may also 
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occur that S vanishes identically, in which case the original equation (3.12| is a conservation law, 
but the expanded form is nonconservative. 

More precisely, within the above setting, we now take into account geometric effects in the 
relativistic Burgers equation (2.11|. By analogy with the first part of this section and in view of 
the equation (3.7), we set x° = ct and x l = r, with c = 1/e and propose the following model: 

/ £ H = ^(-l + V / l + e 2 ^ 2 ), 



(3.13) d t {r 2 w)+d r (r(r~2m)f e (w)^ 



which we refer to as the relativistic Burgers equation on Schwarzschild spacetime. The 

choice of flux and source-term in (3.12) made here allows us to establish a clear connection with 
the model (3.7) derived in the previous subsection from the physicaly sound assumption of Lorentz 
invariance. 

Next, let us seek for steady solutions w s = w s (r), satisfying 



d r (r{r - 2m) f e (w) 







and recalling that f e is strictly convex, we can define its inverse /~ by selecting the "positive 
branch", and we find w s (r) — /~ x f r ( r j^m) ) wnere K i s constant, so that all steady solution to 
(3.13) are given by 



(3.14) 



w a {r) = ± 



K 



r(r — 2m) 



^r(r — 2m) 
J K ' 



3.4. Summary. We have thus derived two Burgers-type models that appear to be of par- 
ticular interest in the context of relativistic fluid dynamics on a curved spacetime. On one hand, 
we have proposed the Burgers equation on Schwarzschild spacetime (3.7), which depends on two 
parameters m, c and reduces to the standard Burgers equation in the limit of zero mass m and 
infinite light speed c. On the other hand, motivated by our derivation of a Lorentz invariant 
model described in the previous section, we have proposed the relativistic Burgers equation with 
geometric effects (3.13), which depends on the light speed 1/e and a geometric weight function UJ. 
Both of these models will be studied in the numerical section below. 

3.5. Alternative approaches. 

Timelike flux vector. Observe that, alternatively, we can normalize T a {v) to be unit time- 
like, g aP T a (v)T f3 (v) = -1, so that T x {v) = ± + e- 2 T°{v) 2 . So, after normalization, we 
find 



(3.15) 



d T° + diV-l + e- 2 ^ ) 2 = 0. 



Derivation from a general covariant equation on Schwarzschild spacetime. We 

return to the analysis in Section |2.5| but now work with the Schwarzschild metric, that is, 



g = -e~ 2 (l - 2m/ r) dt 2 + (1 - 2m/r)" 1 dr 2 
with (a; , a; 1 ) = (t,r). We require that the flux vector field (T°(v), T 1 (v)) is of unit norm, i.e. 



2 (1 - 2m/r) (T») 2 + (1 - 2m/r)" 1 (T 1 ^)) 2 



which yields 



1\2 



(T 1 ) 



2m \ 




Finally, in term of w :— e 2 T°, the associated equation can be written as 



(3.16) 



8qW 



i - 



2m 



1/2 



1 + e 2 1 - 



2m 



1/2 N 



= 0, 
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which, in the limit m — > 0, reduces to the model (2.11) derived in flat space. Then, we can 
introduce z 2 := (1 — 2m/r)w 2 and arrive at 

,3,7) - " /2 ,) + i a ( - 1 + (i - " 2 yrr^) = o. 

4. Well balanced finite volume approximations with geometric effects. 

4.1. Geometric formulation of finite volume schemes. We follow Amorim, LeFloch, 
and Okutmustur [1] and formulate the finite volume scheme for a hyperbolic balance law posed 
on a curved spacetime (M, w), i.e. 

(4.1) divUT(v)) = i (d Q (ojT°(v,x)) + ^(^'(v))) =5(«,x), 

where w : M — > K is the unknown function and T(u) = (T a (v)) and S'(w) : M — > K are pre- 
scribed flux and source fields, respectively, while lo = u>(x) is a positive weight function which 
is identified with the volume form w dx prescribed in some (fixed) global chart of coordinates. 
More precisely, we work on a (1 + 1) -dimensional spacetime and assume the global hyperbolicity 
condition mentioned in the introduction. 

The finite volume methodology applies to general triangulations T = Uifer' 1 ^ °^ ^ ne man i~ 
fold M which are made of spacetime elements K C M satisfying the following structural conditions: 

• The boundary dK of an element K is piecewise smooth dK — Uec9-R" e ano - contains 
exactly two spacelike faces (in the sense (1.3)), denoted by ej- and e^, and "timelike" 
elements 

e°€d°K :=dK\{e+,e K }. 

• The intersection K D K' of two distinct elements K, K' is a common face of K, K' . 

It will be convenient to denote by \K\, |et|, |e^-|, |e°| represent the (one-dimensional or two- 
dimensional) measures of K, e^,e^-,e°, respectively. 

We define the finite volume approximations by formally averaging the balanced law (4.1 ) over 
each element K £ T h of the triangulation. By integrating in space and time, we can write 



div u (T(v)) oj = I S(v)u, 
k J K 

where uo = w(-, ■) is here regarded as a two-form in M. With the Stokes formula, we find 
T°(v)oj(n ep -)= [ T»u;(n e -,.)- £ f T» oj(n e o, ■) + f S(v) u. 

Je K " eOed°K Je " JK 

Given averages and o computed along e^- and e° £ d°K, we need to compute the average 
v\ along e\. To this end, we introduce 



T°(v)u)(n e +,-) - \e K \w e - T - (v K ), / T 1 (v) u(n e o , ■) ~ \e°\w e o Q Kte o(v K ,v K ), 

and 

/ S(v)uj ~ \K\uj k S K , 

JK 

in which, for each K and e° £ d°K we have selected a (Lipschitz continuous) numerical flux 
Qk.c : K 2 — > R satisfying natural properties of consistency, conservation, and monotonicity. In 
turn, the finite volume method of interest takes the form 

(4.2) \e^\u} e +T e +(v^) = \e^\uJ e -T e -(v^)- l e °l w e° <3A-,e»(%.\„) + \ K \ ^kS K - 

e°ed°K 

For the sake of stability of this scheme, a standard CFL condition is assumed. See [I] for further 
details. 

H 



4.2. Finite volume schemes in coordinates. From now on, motivated by the models 
discussed in the previous two sections, we assume that the spacetime is described in coordinates 
(t, r) such that the weight oj = w(r) as well as the flux and source terms depend only upon the 
spatial variable. We divide into equally spaced cells Ij — [r j_ 1/2, Tj +1/2] of size Ar, centered at 
Tj, that is, fj+1/2 = r j-i/2 + Ar, satisfying Tj-\i 2 = jAr, rj = (j + 1/2) Ar. Moreover we denote 
by At the constant time length and we set t n = nAt. 

As just explained, the finite volume method is based on averaging the balance law (4.1) over 
each grid cell [t n ,t n +i] x Ij, which yields 

f diy u (T(v))cJdx°dx 1 = [ S(v)uJdx°dx 1 

J/jX[t„,t„ + i] Jl :i x[t n ,t n+1 ] 

and rearranging the terms gives 

j S(v)ujdx°dx 1 = f 3+1/2 (T°(t n+1 ,r)-T°(t n ,r))uJdx 1 

J/ 3 x[t„,t„ +1 ] •'r j _i/2 

/tn+l 
(aJ(r i+1/2 )T 1 (<,r J+1/2 ) -uj{r j _ 1/2 )T 1 (t,r j _ 1/2 ))dx . 

We recall that uj is independent of the time variable, so we can write 

■Tj+i/2 m+1/2 r 

T^tn+^^Lddx 1 = T^tn^^Ujdx 1 + S^Zudx^dx 1 

Tj-1/2 J IjX[t n ,t n + 1 ] 

-fwj +1 /2^ T x (t, r j+1/2 )dx° -uj j _ 1/2 J^ T 1 (t,r j - 1 / 2 )dx ^. 
We next introduce the following approximations of numerical flux functions and the source term 
4t T\t,r J±1/2 )dx° ~ Q? ±1/2 , j- r +1/2 T*{t n ,r)udx l ~ T?^ 

and 

-j — 7— / S(v)uJdx°dx 1 ~ S™ 5J,-. 

Hence, based on these approximations, the finite volume approximation takes the form 



^{^] + l/2 Qj + 1/2 - Uj-l/2 Qj-l/2) + AtSj Wj, 



where A := At/Ar. Since Tj = T(v™) is invertible, after dividing the above equation by Uj, we 
arrive at the scheme 



(4-3) v] +1 = T 1 (t(^) - - (m j+1/2 Q] +1/2 oJ^x/2 QU/2) + At S{v?)j . 

4.3. Well balanced scheme for hyperbolic equations on Schwarzschild spacetime. 

We can think of a general solution to ( 4.1 ) as being, locally in space and time, a small perturbation 
of a steady solution, that is, in a regime where the source term balances the flux gradient within 
a certain moving frame. Asymptotically in time, one expects the behavior to be steady, since the 
geometry is static, as is the case for the models in the previous section. We focus attention on 
defining a well balanced scheme specifically in the case 

ui(r) = r(r — 2m), 
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for some mass m > 0, and we present the scheme in the case of the Burgers equation on 
Schwarzschild spacetime, that is, (3.7). The discrete version of (3.7| now reads 

(4.4) 



7ri n + 1 _ 7p n 

3 3 



At 

Ar 



( UJ 3 + l/2Q] + l/2 - Wj-l/iQj-l/i) +AtSj, 



but we now take into account the mass 



where Ar still denotes the mesh size Ar = r J+1 / 2 — rj-i/21 
parameter in order to define the mesh points, that is, 

fj-i/2 = 2m + jAr, rj = 2m + (j + l/2)Ar, r j+1/2 = 2m + (j + l)Ar 

and the averaged weights are Wj±i/2 = i~j±i/2( r j±i/2 — 2m). Moreover, T™ and Sj are approxima- 
tions corresponding to 



T 3 = 



1 



AV Jr., 



r 3 + l/2 _ 1 

TdV~ gi S 3 = 



3-1/2 



r 3 + l/2 



AV Jr j _ 1 



SdVg, 



in which g is the induced metric on spacelike slices of the foliation, with g = —c 2 (l — dt 2 + g 
and dVg is the natural volume form induced on a spacelike slice, so 



dV~ g := r 2 



V r — 2m J 



1/2 



dr, 



AV := r 2 



r — 2m 



1/2 



Ar. 



It is convenient to introduce Vl := luv 2 /2. To treat the source term, we integrate the equation 



(3.8) over a spacetime cell, and obtain 
rr+1/2- 



' g Q,dr = 



r-l/2+ 
rr+1/2- 



r+1/2- 



r—l/2+ 



[rv 2 — mc 2 ) \fgdr — AVSj, 



which gives Sj = J^_ 1 j 2+ yy^dr. By substitution and integration by parts, we arrive at the 
expression of the source term: 



S } = 



-L(r 7 / 2 (r-2m) 1 ' 2 v 2 /2) 



-+1/2- 



1 



-1/2+ AV J T -\/2+ 



r+1/2- r ( r _ 2m ) 



r r(r — 2m) 



dVs 



In the numerical experiements, we will use Gaussian quadrature to evaluate this integral. 
Furthermore, the left- and right-hand numerical flux terms are given by 



Qj+1/2 - Q(Tj+i/2-,T j+ i/ 



2+ h 



0,-1/2 = Q(T 



j-l/2-,Tj-l/2+), 



where Tj + i/ 2 ± = r j+i/2 v j+i/2± with reconstructed velocities 



v j+1/2+ = ±Jc 2 -K 2 [l 



2m 

r j+ i/2 



v^ 1/2+ =±Jc 2 -K 2 [l- 



r j-l/2 



U i + l/2- 



±Ac 2 -K 2 \\ 



fl-— ) 

r 3 + l/2 J 



2m 

r 3-l/2 



and 



9 —9 

c 2 - v 2 +1 



K 



1 —9 
" - V 3 



Ki 



9 —9 

c 2 - v 2 _ x 



This reconstruction is motivated by the steady solutions derived in (3.10) which allow us to 
compute intermediate states at which the numerical flux is evaluated. Importantly, in the special 
case that K \ , K 2 , K3 coincide, the proposed scheme recovers an exact steady solution. 
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4.4. Second-order accurate version. To increase the order of accuracy of the above 
scheme, we rely on the technique introduced by Nessyahu and Tadmor central scheme [5]. Our 
second-order version is thus a predictor-corrector scheme, defined by 



- \ K + i + -U) + - *W - 2 ( T >m /2 ) - 

where A := At/ Ar denotes the mesh ratio and 8qT^ j is the approximate first-order derivative 
of the flux T^ j = The time step must satisfy the stability condition A max^ |^)T^| < |, 

where Ido? 1 * -| denotes the spectral radius of the Jacobian matrix. In order to get rid of spurious 
oscillations, a slope limiter must be applied to the values v'j and doT^j and, specifically, 

Vj = Minmod (wj+i — Vj, Vj — 
doT^j = Minmod (T/ +1 - T/, T] - T}_ x ), 
Minmod(a, b) = sgn(a) min(|a|, if ab > 0; if ab < 0. 

5. Numerical experiments. 

5.1. Comparison between several schemes and models. Numerical experiments are 
now presented for both models derived in this paper. After normalization (c = 1/e = 1), wc thus 
consider the geometric Burgers equation (I): 

(5.1) d t (r 2 w) + d r (r(r- 2m) (- 1 + ^l + w 2 )^j = 0, 
and the geometric Burgers equation (II): 

(5.2) d t (r 2 v) + d r (r(r - 2m) V -] = rv 2 - m, 



with obvious notation. The first equation has a conservative form, while the second one is non- 
conservative. We work within the interval r € (2m, R), where R is a fixed upper bound for the 
spatial variable. We impose no-influx boundary conditions at r = 2m and r = R by solving a 
Riemann problem between the boundary data determined by the initial data and the current nu- 
merical value at the boundary. In other words, we use the Godunov flux at the boundary. For all 
numerical experiments, we use 2m = 0.1 and R = 1.0 and the value of CFL number is taken to be 
0.9. We begin by illustrating the importance of the well balanced property by comparing together 
the numerical solutions obtained with the proposed scheme and with a more naive discretization of 
the right-hand side of the geometric Burgers equation (II). For the model (I), we used Ar = 0.02 
and the endpoint velocity v(R) — 0.03. The numerical solution is plotted on Figure 5.1.1. In this 
Figure, in (b), (c), we plot the numerical solutions given by the three schemes. The first one is the 
first-order Lax-Friedrichs scheme (which is indicated by +); the second one is the second-order 
version of Lax-Friedrichs scheme (which is indicated by line); the third one is the well balanced 
second-order version of Lax-Friedrichs scheme (which is indicated by — ). In Figures (d), (e), we 



compare the numerical flux r(r — 2m) f e (w) (cf. (5.1 1), by the three schemes. The proposed well 



balanced scheme appears to be efficient and robust and, in particular, this scheme preserves steady 
solutions. Next, for the model (II), we used Ar = 0.005 and the velocity v(R) — 0.3, which led 
us to Figure 5.1.2. In Figures (b) and (c), we compare the numerical solutions based on the three 
schemes. The first one is a standard first-order Lax-Friedrichs scheme (which is indicated by 
+); the second one is a standard second-order Lax-Friedrichs scheme (and is plotted with dots); 
the third one is a well balanced second-order Lax-Friedrichs scheme (and is plotted with — ). In 



Figures (d) and (e), we compare the numerical values of K 2 (defined in ( |3.10[ )) based on the three 
schemes. Again, the scheme is found to be efficient and preserve steady solutions. 
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(a) 




5.2. Comparison of several schemes for a single shock. We study next weak solutions 
containing a single shock separating two steady solutions at the location r = 0.5: one steady 
solution corresponds to the data Vi(R) =0.1 while the other steady solution corresponds to 
V2{R) — 0.9. The initial data are constructed so that the initial discontinuity propagates as a single 
shock (without decomposing itself during the evolution), which interacts with the background 
geometry. We consider the second model (5.2) with Ar = 0.05 . The numerical solution is 
presented in Figure 5.2. We compare the propagation of a single shock by three schemes. The 
first one is a standard first-order Lax-Friedrichs scheme (plotted with a continuous line) , and the 
second one is a standard second-order Lax-Friedrichs scheme (ploted with the symbol +). Finally, 
the third scheme is the proposed well balanced second-order Lax-Friedrichs scheme (plotted with 
the symbol x). One observes the efficiency and robustness of the well balanced scheme; the 
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(a) 




numerical solution based on this scheme is monotone and stable under perturbations. 

5.3. Application to late-time asymptotics perturbed steady solutions. We study 
here the nonlinear stability of a steady solution superimposed with a compactly supported initial 



perturbation. In Figure 5.3.1, we still use (5.1 ) and take each subinterval Ar = 0.03 and a velocity 
at the right-hand point v(R) = 0.03. We observe that the initial perturbation evolves toward a 
so-called iV-wave, which is typical of nonlinear hyperbolic systems, and that this perturbation 
moves in the right direction away from the singularity. It eventually goes away at the boundary 
r = R of the computational domain. We observe the convergence of the solution toward the 
steady solution using the well balanced second-order Lax-Friedrichs scheme. One observes that 
the numerical solution, using the well balanced second-order Lax-Friedrichs scheme, is stable. 
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(a) 



(b) 





(c) (d) 
Figure 5.2. The numerical solutions given by the three schemes for the second model (II). 



Next, in Figure 5.3.2, we treated the second model (5.2), with Ar = 0.001 and v(R) = 0.1, and 
similarly to the first case, one observes that the initial perturbation evolves toward a so-called 
iV-wave, and that this perturbation moves in the right direction away from the singularity. Again, 
the solution converges toward the steady solution. 

5.4. Application to late-time asymptotics perturbed one shock solution. Finally, 
we take the case of a solution containing a single steady shock separating two initially steady 
solutions for the equation (3.14) at the location r = 0.5. We consider here the first model (5.1) 
with Ar = 0.05 and v(R) = 0.09. We study here the nonlinear stability of this steady solution and 
include initial perturbation on a steady solution. The numerical results are plotted in Figure 5.4. 
We observe the convergence of the solution toward the steady solution (single steady shock) using 
the well balanced second-order Lax-Friedrichs scheme. 

6. Conclusions. In the first part of this paper, we considered a general class of hyperbolic 
balance laws posed on a curved spacetimes and we derived flux vector fields T = T{v) which, in 
a relativistic context, generalize the classical Burgers flux. On one hand, we required that (1.1) 



satisfies the Lorentz invariance property shared by the Euler equations of relativistic compressible 
fluids, and we identified a balance law, which is unique up to normalization. Interestingly enough, 
the standard Burgers equation was recovered by taking the singular limit of infinite light speed. 
On the other hand, an alternative model of interest was derived directly from the relativistic 
Euler equations, by assuming that the pressure term vanishes. This latter model, in the case 
of a flat background geometry, turned out to be simply the classical Burgers equation. Both 
models are referred to relativistic Burgers equations on curved spacetimes. They exhibit many of 
the mathematical properties (hyperbolicity, genuine nonlincarity, steady solutions) and challenges 
(shock wave, late-time asymptotics) encountered with the Euler system of relativistic compressible 
fluids. In the second part of this paper, we investigated the numerical discretization of entropy 
solutions to the above two models. We introduced a finite volume scheme for ( |1.1[ ) defined on 
curved spacetimes, especially on a background Schwarzschild spacctime. The proposed scheme is 
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(a) 



(b) 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(c) 



(d) 



Figure 5.3.1. Model (I). The asymptotic behaviors of long time solutions for equation (5.1| using 

the well balanced second-order scheme. 



fully consistent with the divergence form of the equations and, therefore, applies to weak solutions 
containing shock waves. More importantly, our scheme is well balanced in the sense that it 
preserves steady solutions, and extensive numerical experiments demonstrated the convergence of 
the proposed scheme, and its relevance for computing entropy solutions on a curved background. 
The generalization to the Euler equations posed on Schwarzschild spacetime are presented in the 
follow-up work [BJ. 
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